- One factor ANOVA
- Git and GitHub
- Means tests in ANOVA
- Experimental Design
- Power analyses
- Begin Multi-factor ANOVA
10/30/2018
git clone https://github.com/wcresko/UO_ABS.git
git status git merge origin/master
dplyr functions.RNAseq_Data <- read.table("RNAseq.tsv", header=T, sep='')
x <- RNAseq_Data$categorical_var
y <- RNAseq_Data$continuous_var1
z <- RNAseq_Data$continuous_var2
contrasts(x) <- cbind(c(0, 1, 0, -1), c(2, -1, 0, -1), c(-1, -1, 3, -1))
round(crossprod(contrasts(x)), 2)
rnaseq_data_list <- list(x = list(‘xxx vs. xxx’ = 1, ‘xxx vs. xxx’ = 2, ‘xxx vs. xxx’ = 3))
RNAseq_aov_fixed <- aov(y ~ x) plot(RNAseq_aov_fixed) boxplot(y ~ x) summary(RNAseq_aov_fixed, split = rnaseq_data_list)
perc <- read.table('perchlorate_data.tsv', header=T, sep='\t')
x <- perc$Perchlorate_Level
y <- log10(perc$T4_Hormone_Level)
MyANOVA <- aov(y ~ x)
summary (MyANOVA)
boxplot(y ~ x)
install.packages("multcomp")
library(multcomp)
summary(glht(MyANOVA, linfct = mcp(x = "Tukey")))
Survival of climbers of Mount Everest is higher for individuals taking supplemental oxygen than those who don’t.
Why?
The goal of experimental design is to eliminate bias and to reduce sampling error when estimating and testing effects of one variable on another.
\[ Power \propto \frac{(ES)(\alpha)(\sqrt n)}{\sigma}\]
Power is proportional to the combination of these parameters
perc <- read.table('perchlorate_data.tsv', header=T, sep='\t')
x <- perc$Strain
y <- log10(perc$T4_Hormone_Level)
MyANOVA <- aov(y ~ x)
summary (MyANOVA)
boxplot(y ~ x)
Based on your results, calculate the power for your ANOVA.
pwr.t2n.test(n1=xxx, n2=xxx, d=xxx, sig.level=.05, power=NULL)
Check out the functions in the ‘pwr’ library (Unnecessary in this case, but could use ANOVA version):
pwr.anova.test(k=2, n=190, f=0.25, sig.level=.05, power=NULL)
andrew_data <- read.table('andrew.tsv', header=T, sep=‘\t')
head(andrew_data)
andrew_data$TREAT2 <- factor(c(rep(“low”,40),rep(“high”,40))
The nested factor is PATCH - also need to turn this into a factor
andrew_data$PATCH <- factor(andrew_data$PATCH)
In this case, our response variable is ALGAE Look at the distribution of ALGAE for the two levels of TREAT2 using boxplots based on the patch means, which are the replicates in this case.
andrew.agg <- with(andrew_data, aggregate(data.frame(ALGAE), by = list(TREAT2=TREAT2, PATCH=PATCH), mean) library(nlme) andrew.agg <- gsummary(andrew_data, groups=andrew_data$PATCH) boxplot(ALGAE ~ TREAT2, andrew.agg)
Run the nested ANOVA:
nested.aov <- aov(ALGAE ~ TREAT2 + Error(PATCH), data=andrew_data) summary(nested.aov)
Do we detect an effect of TREAT2 (high vs low sea urchin density)?
Estimate variance components to assess relative contributions of the random factors
library(nlme) VarCorr(lme(ALGAE ~ 1, random = ~1 | TREAT2/PATCH, andrew_data))
rnadata <- read.table('RNAseq.tsv', header=T, sep='')
head(rnadata)
continuous response variable and two main effect categorical variables
gene <- rnadata$Gene80 microbiota <- rnadata$Microbiota genotype <- rnadata$Genotype boxplot(gene ~ microbiota) boxplot(gene ~ genotype) boxplot(gene ~ microbiota*genotype)
rna_aov <- aov(gene ~ microbiota + genotype + microbiota:genotype) rna_aov <- aov(gene ~ microbiota*genotype)
Examine the fitted model diagnostics and the ANOVA results table
plot(rna_aov) summary(rna_aov) anova(rna_aov)
rnadata <- read.table('RNAseq.tsv', header=T, sep='')
head(rnadata)
variables excluding first 5 and last 5 observations
gene <- rnadata$Gene80[6:75] microbiota <- rnadata$Microbiota[6:75] genotype <- rnadata$Genotype[6:75] boxplot(gene ~ microbiota) boxplot(gene ~ genotype) boxplot(gene ~ microbiota*genotype)
Estimate the variance components using Restricted Maximum Likelihood (REML)
library(lme4) lmer(gene ~ 1 + (1 | microbiota) + (1 | genotype) + (1 | microbiota:genotype))
Based on the REML sd estimates, what are the relative contributions of the factors to total variance in gene expression?
longevity_data <- read.table('longevity.csv', header=T, sep=',')
head(longevity_data)
Variables
long <- longevity_data$LONGEVITY treat <- longevity_data$TREATMENT thorax <- longevity_data$THORAX
boxplot(long ~ treat) plot(long ~ thorax)
plot(aov(long ~ thorax + treat ), which = 1)
plot(aov(log10(long) ~ thorax + treat ), which = 1)
library(lattice)
print(xyplot(log10(long) ~ thorax | treat, type = c("r", "p")))
anova(aov(log10(long) ~ thorax*treat))
anova(aov(thorax ~ treat))